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We study the sound velocity in cubic and non-cubic three-dimensional optical lattices. We 
show how the van Hove singularity of the free Fermi gas is smoothened by interactions and 
eventually vanishes when interactions are strong enough. For non-cubic lattices, we show 
that the speed of sound (Bogoliubov- Anderson phonon) shows clear signatures of dimensional 
crossover both in the ID and 2D limits. 
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I. INTRODUCTION 

The experimental ability to tune and control the quantum properties of matter has increased 
dramatically during the past decade. The evaporative cooling of bosons enabled the creation of 
Bose-Einstein condensates in dilute atomic gases . Due to statistics fermionic atoms proved 

harder to cool, but the favorable collisional properties of certain two-component mixtures made 
the efficient cooling of fermions into quantum degeneracy possible P, Q|. 

If two different kinds of fermions are mixed, they can interact via s-wave interactions. If this 
interaction is attractive, fermions can pair and form a superfluid as predicted by the Bardeen- 
Cooper-Schrieffer (BCS) theory [7j. However, the transition temperature of the superfluid tran- 
sition is typically much smaller than degeneracy temperature of the fermions and therefore hard 
to reach experimentally. This difficulty was recently overcome by using the so-called Feshbach 
resonance to change the strength of the interaction between two atoms. 

Near a Feshbach resonance an external magnetic field 8] is used to tune the energy difference 
between the molecular bound state and the two-atom continuum. Changes in this energy difference 
give rise to a resonant behavior of the interaction strength between two atoms. In particular, the 
interaction strength can be varied from a small and negative value (BCS-side), into a very large 
value close to resonance, and finally into a small and positive value (BEC-side) at the other side 
of the resonance. In the Bose-Einstein condensate side of the resonance the system forms diatomic 



molecules 



111 ] and, at sufficiently low temperature, a condensate of molecules |l2l.ll3l.ll4l.ll5|. 
Fermion pairing is possible on the BCS-side of the resonance. Strong indications of such pairing 
were experimentally observed employing magnetic field sweeps across the resonance |lfit Il7| and 
by monitoring the behavior of collective modes 
measured using radio-frequency spectroscopy |2£ 



Q|. Soon alter the pairing gap was directly 
Finally, the smoking gun of superfluidity, a 



vortex lattice, in a fermionic system was very recently observed 

New tools to manipulate the atomic clouds have been developed. In particular, optical lattices 
have proved to be an especially useful tool. In the context of degenerate quantum gases, optical 
lattices 0, Q have been used to explore, for example, non-classical states, phase co- 
herence, condensate dynamics, as well as matter wave interference. Since both the lattice depth 
and the lattice geometry can be easily changed, there is the possibility of investigating experi- 
mentally anisotropic quantum systems as well as quantum phase transitions from, for example, a 



Q. 



superfluid into an insulator as the tunneling strength between nearest neighbors is reduced 
Also, sufficiently deep lattices in certain directions can be used to control the dimensionality of the 
system. 

The ability to trap atoms into optical lattices does not depend on which statistics, fermionic or 
bosonic, the atoms obey. Indeed, recently degenerate fermions have been studied experimentally 



in one-dimensional 



231 as well as in three-dimensional 



3,0 



optical lattices. There has also 



been rapid progress on the theoretical front. Among other 



32 



things, there are several studies on 



33 



34, 



35, 



36, 



371 ] and on fermion 



fermions and boson-fermion mixtures in optical lattices 
dynamics in optical lattices In addition, more exotic fermionic systems with several 

flavors have been discussed [4l(. 

In this paper we calculate the density response of superfluid fermions in a three-dimensional 
(3D) optical lattice, to elucidate how to investigate dimensional crossover effects. Therefore we 
calculate the above signatures also for non-cubic optical lattices with varying lattice depths in 
different directions. We show that the speed of sound (Bogoliubov- Anderson phonon) as a function 
of filling fraction behaves qualitatively differently when the dimensionality of the lattice is changed. 

While we work consistently within the framework of the BCS theory, we expect our results and 
discussions to be instructive also somewhat closer to the Feshbach resonance where interactions 
are stronger. Naturally it is to be understood, that if the scattering length characterizing the 



interaction strengt 



will change 



42j,|43 



r becomes comparable to the lattice spacing, the continuum Feshbach physics 
If we furthermore have about one atom per lattice site, the simple BCS 



description is expected to fail in this case. 



This paper is organized as follows. In sec. [H] we review the standard BCS-theory. In section 
II I II we discuss the density response and how it is calculated. In section HVl we present and analyze 
our results on the speed of first sound in cubic lattices and in section study the dimensional 
crossover in non-cubic lattices. Conclusions are presented in sec. IVII 



II. BCS-THEORY 



The single-band Hubbard Hamiltonian in the context of ultracold gases can be derived from 
the full microscopic Hamiltonian j^|. To achieve this, one must expand the field operators of 
the microscopic Hamiltonian in the basis of the well localized Wannier functions, assume that 
the optical lattice is sufficiently deep, and that the gas is sufficiently dilute. Then it is possible 
neglect all but the nearest neighbour interactions. Furthermore, assuming that the temperature 
and interactions are low enough so that there are no excitations to the excited oscillator levels of 
the individual lattice wells, the Hamiltonian becomes, for fermions, the single-band Fermi-Hubbard 
Hamiltonian: 

H - llN = - fl ^ (4rAn + c\ n Cln) + U ^ c\ n c\ n C ln C^ n 
n n 

( \/ N 

~ Ux Yl +J V Yl + Jz Yl I (4™ a Tn + c\mhn) , 
\ (n,m) x {n,m) y (n,m) z J 

Here the operator ctn creates an atom with spin a at the lattice site n and c an annihilates it. 
The bracket (n, Tfi) & denotes all the nearest neighbour pairs in the ^-direction, ^jl is the chemical 
potential, U is the interaction strength of the atoms and J x , J y , and J z are tunneling strengths 
in their respective directions. We assume that there is an equal amount of spin up and spin down 
atoms, so the chemical potentials for both states are equal. In order to allow the lattice to be 
non-cubic, separate values for the three tunneling strengths are needed. 

The parameters J and U can, in principle, be computed numerically using the lowest Wannier 
function, but if one approximates the lattice well by a harmonic potential they can be determined 
analytically. The interaction parameter is given by 

8 a ( 2s x s y s z E 3 R h 2 y 



while the hopping parameter is given by 



[4^ 



./,-^-^(('^i-^)-i Si (l + e -v^)). 



Here A is the wavelength of the lattice, o is the scattering length, = fr? (2tt / A) 2 / '2m is the 
recoil energy of the lattice and s x , s y and s z are the lattice heights, in recoil energies, in different 
directions. 

The energy separation of the states of individual lattice wells is 

where s is the lattice height. Prom this and formula Eq. © it is possible to estimate the limits 
of validity for the single band Hubbard model, i.e. when thermal excitations or interactions are 
strong enough for the second band to be populated. For 6 Li with A = 1030 nm and s = 2.5 these 
limits are approximately 3 [iK for the temperature and —8600 ao for the scattering length. This 
upper limit for the scattering length is so high that typically the on-site interaction approximation 
fails before reaching this limit. 

The Hamiltonian in momentum space is 

H = Yl ( ek ~ (4iAk + 4k £ ik) + ^ H 4k'+k+ q 4k'-k- q 2 !k'+kc T k'-k, (5) 

k k,q,k' 

where M is the number of lattice sites and the lattice dispersion is given by 

e k = 2J X (1 - cos(k x d)) + 2J y (l - cos(k y d)) + 2J Z (1 - cos(M))- (6) 

Note that the dispersion is even, so = e_k. The Heisenberg equations of motion for the operators 
c-f_k and 4k are 

-ih— c T _ k = [#,c T _ k ] 

a ( ? ) 

-^4k = t^4k]- 

As it stands, the Hamiltonian in Eq. @ also contains pairs with net momentum. According to 
the BCS approximation, such terms can be dropped. This is formally achieved by setting k' = 
in Eq. Q- With this simplification, the commutators are 

[H, c T - k ] = - (ek - M)q-k - -jj^cjk Yl £ !-k-qCTk+q 



[^4k] = ( £ k " M)4k "fE (4k+q4-k-q) fi T-k- 

q 

These equations are then linearized by replacing the sum £V cj.-k-qC|k+q with its expectation 
value. When the chemical potentials for different components are the same, this value is always 



(C|-k-qC T k+q) = A 

q (9) 
-V/c* gt \ =A *. 

M ^ \ Tk+q |-k-q/ 



independent of position in equilibrium. 

U_ 

~M 

U_ 

M 
q 

Because the overall phase does not play any role in the mean-field theory, it is possible to select 
the order parameter A to be real. The equations of motion are thus reduced to 

These equations are decoupled with the standard Bogoliubov-transformation, i.e. with a transfor- 
mation that diagonalizes the Hamiltonian 



C|-k \ / u-k «-k 



4k I \ _, 'k "k 




(11) 



For this transformation to be canonical, the new operators must obey the fermionic anticommuta- 



tion rules. This leads to the requirement that u? + 1>? = 1. The transformation that decouples the 



equations is given by 



2 1 /, . e k - M 

K = - 1 1 + 



2 V #k 
2 V #k 



4 = ^(1-^) (12) 



nkUk= 2ib 

where -E^ is the quasiparticle dispersion, i.e. 



E k = V(e k -M) 2 + A2. (13) 
Inserting the transformed operators back to Eq. gives the so called gap equation: 

A = H (c|-kC T k> =-TjJ2 \ (- w -k7{k + M -k7i-k) • (^k7tk + ^k7i_k) )• ( U ) 
k k ^ ' 

Because the quasiparticle operators j describe fermions, they obey Fermi statistics. Therefore, in 
thermal equilibrium, 

(7| k 7 T k> = f(E k ) 
<7|-k7l_ k > = 1 " f(E-k) 

(nr&U) = 

(7i-k7?k> = 0, 



and 

A=4s4 (1 - 2/(£k)) ' (16) 

where / is the Fermi distribution. On the other hand, the number of particles is equal to 

N = E (4kCik + 4 k c?k) (17) 

k 

and using Eq. (|12|) and Eq. (|15|) this becomes the number equation, 

N = 2 £ K/(S k ) + - /(£ k ))) . (18) 

k 

Instead of the total number of particles, it is useful to deal with the filling fraction, i.e. the total 
number of particles divided by the number of lattice sites, N/M. In our notation this ratio can 
have values between and 2. The former corresponds to no particles, whereas the latter describes 
a full lattice with two atoms of opposite spins at each site. 



III. DENSITY RESPONSE 



Density response describes how the total density of the system changes as a result of a (sharp) 
change in the external potential. In the linear response regime, in a homogenous systema and 
under equilibrium conditions, it is possible to write 

6pQt,u,)=xQL,u)SVQL,u>), (19) 

where 5V(\c, to) is the change in the external potential and x(k, w) is the density response. 

To calculate the response function x(k, u>), we use the method of the generalized random phase 
approximation (GRPA) following the notation used by Cote and Griffin 49]. Let us 

briefly summarize the final results of their derivations. First, define two matrices, A and B: 



*<k)-l 3k "" k ! k |. Pi) 



Then define a 4 x 4-matrix, L°: 



^1111 


^1121 


^1211 


L 


^1112 


£ll22 


^1212 


L 


^2111 


-^2121 


^2211 


L 



2221 



(22) 



\-^2112 -^2122 -^2212 -^2222/ 



with matrix elements given by 

l ij (q + q / )Ifcz(q / ) 



Eq> - -Eq+q' + H(w + iS) 



(/-/') + E ^ +q !t (q ?^ (/ 



^q' 



q+q' 



/') 



^(q + qQBfcKqQ , , _ , g^q + qQ^fcKqQ , , 



(23) 



E a / + .E, 



q+q' 



h(u) + i5) 



Eq 1 + -Eq+q' + ^(uJ + 15) 



1) 

rence factor 



where / is the Fermi function at E„i, f is the Fermi function at -E q '_ q , and the converg 
S is put to after the calculation. 

Due to symmetries, only six of the 16 elements in L° are actually independent |49|. These 
elements can be denoted as 



a = 


£1111 


= £2222 






b = 


£1212 


= —-^1221 


= — £2112 


= £2121 


c = 


£1112 


= —-^1121 


= L1222 = 


— £2122 


c = 


£1211 


= — £2111 


= £2212 = 


— £2221 


d = 


£1122 








d = 


£2211 









(24) 



The number of independent components can be reduced further by assuming weak coupling, as 
is done in j^jj. In this limit, c = — c and d = d. However, in order to probe also more strongly 
interacting systems, we do not make this assumption. We then define a 4 x 1 vector L, 

(a-b\ 
2c 
2c 
\a-bj 

and solve the linear algebra problem 



L 



(25) 



(1 + UL°) x = L, 



(26) 



FIG. 1: The Anderson-Bogoliubov phonon in an cubic lattice. The x-axis is the x-component of the 
momentum vector in the units of inverse lattice spacing, the y-axis is energy in the units of the pairing gap 
divided by H, and the z-axis shows the dynamical structure factor from Eq. Q28JI. in arbitrary units. 

where U is the interaction strength in the Hamiltonian. Finally, from the solution x we deduce the 
density response function from x: 



We consider the dynamic structure factor instead of the e xp 



icit density response, because the 



former is measurable using, for example, Bragg spectroscopy [42|. The dynamic structure factor is 
given by 

S(k,Lo)=--lm X (k,uj). (28) 

7T 

Formally, the convergence factor 5 denned in Eq. Q23|) is put to after the calculations. However, 
for the sake of illustration, we use a small finite value for 5 in the figures. This gives rise to a finite 
linewidth in the density response. 

Because the experimentally relevant lattice sizes are relatively small, only in the order of 100 
sites per dimension, it is possible to explicitly calculate, i.e. without approximating by integrals, 
the sums over all lattice sites in the BCS-theory and the density response. For concreteness all our 
results are calculated for 6 Li atoms at zero temperature and with a laser wavelength A = 1030 nm. 
Restriction into T = is accurate when the energy gap is much larger than the temperature. 
Our numerical calculation correctly reproduces the Anderson-Bogoliubov phonon, see Fig. ^ This 
mode is gapless and is expected on general grounds when the continuous U(l) symmetry is broken 
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FIG. 2: The dispersion of the Anderson-Bogoliubov phonon. The x-axis is the x-component of the momen- 
tum vector in the units of inverse lattice spacing and the the y-axis is energy in the units of the pairing gap 
divided by H. 



The density response in (k, uj) space (see Fig. |2J) gives the dispersion of sound in the lattice. The 
long wavelength limit of the dispersion is linear and therefore the speed of sound is independent of 
momentum in that limit. For higher momenta the response saturates to 2A. In general the speed 
of sound in a weakly interacting Fermi gas is (to the leading order) vf/V3, where vf is the Fermi 
velocity of the system. However, since the Fermi surface in a lattice is not a sphere, the Fermi 
velocity is not a uniquely defined number, but depends on the direction. With the definition 



our results agree with the free space result in the low density limit. 

In order to experimentally observe the density response in a finite system one must overcome 
some problems. In particular, the region of k-space where the dispersion relation is linear scales 
with the inverse of the coherence length (or Cooper pair size) i.e. 5k ~ l/£. In the weakly 
interacting system, Cooper pairs become very large and the density response saturates very quickly 
to 2 A. In the lattice the smallest non-zero wavevector has the magnitude of fe™ m = 2ir/N x \, where 
N x is the lattice size in the x-direction. This should be much smaller than the linear region of 
the dispersion. In other words the system size must be much larger than the coherence length. 
This leads to a condition for the minimum system size where the observation of density response 
is possible, which is shown in Fig. [3] 



IV. SPEED OF SOUND IN A CUBIC 3D LATTICE 
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FIG. 3: Pair size compared with the system size. The x-axis is the number of the lattice sites in the re- 
direction and the y-axis is the absolute value of the scattering length, in Bohr radii. For the density response 
to be experimentally observable, the coherence length, i.e. the size of a Cooper pair, must be smaller than 
the system. The linear sound dispersion is clearly visible above the solid line. 



Unlike in free space, in a lattice the speed of sound changes with the density of the atoms in 
a non-trivial way. Fig. 0] shows the speed as a function of the filling fraction with different inter- 
action strengths. The sound velocity is typically on the order of mm/s and is thus experimentally 
measurable. With small scattering lengths and in a cubic lattice, the speed of sound as a function 
of the filling fraction is non-monotonous and exhibits a local minimum around the point where 
Ep = 4J. For non-interacting fermions this minimum is a sharp cusp corresponding to the van 
Hove singularity in the density of states. The cusp appears when the Fermi surface in the direction 
of the sound reaches the divergence in the density of states. The cusp is visible in Fig. 0] where we 
show the isothermal sound velocity 



where P is the pressure and p is the mass density of the gas, calculated for the ideal Fermi gas in 
a finite sized lattice. 

However, interactions smear the Fermi surface by broadening its edge by an amount of the 
pairing gap, A. This smoothens out the cusp caused by the singularity. In fact, as is clear from 
the Fig. |1J the effect will vanish completely with strong enough interactions. 

The finite size effects introduce some modulations, but it is clear that a non-zero gap shifts 
the minimum of the sound velocity into higher filling fractions. The non-zero gap also lowers the 
sound velocity and this effect becomes more pronounced at higher filling fractions. For a weakly 



A. Speed of sound and the van Hove singularity 
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FIG. 4: The speed of sound as a function of the lattice filling fraction, for the non-interacting Fermi gas (top 
curve) and for interacting Fermi gases with different scattering lenghts: (in Bohr radii, from top to bottom) 
—750, —1000, —1200, —1500, and —2000. The speed has a local minimum at a point where Ep — 4 J, but 
this minimum is smoothened out as interactions become stronger. 
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FIG. 5: The Fermi surface of the lattice with different densities, as seen from the z-direction. The axis are 
the x- and y-components of the momentum vector. The surfaces above were plotted with filling fractions 0.2, 
0.44, 0.47, and 0.8, from left to right. The colour coding: white corresponds to the Fermi surface reaching 
the edge of the Brillouin zone in the z-direction, whereas the black area is outside the surface. Note that 
the presence of the excitation gap spreads the Fermi surface, wheras for a non-interacting system, the edges 
of the Brillouin zone are not reached until the filling fraction exceeds 0.42. 



interacting system with a very low filling fraction and an infinite lattice, the speed of sound is given 



by the analytical result 



48 




[/(67r 2 n) x /3 



4vr 2 J T d 2 



(31) 



The location of the van Hove singularity is also the point where the Fermi surface of the system 
becomes disconnected in a non-interacting system. In the interacting case, the Fermi surface is 
smoothened, see Fig. 03 
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FIG. 6: The speed of sound as a function of the scattering length a. Here d is the lattice spacing, the filling 
fraction is 0.5, and the lattice is cubic, with s x , s y , s z = 2.5. 

B. Speed of sound as a function of interaction strength 

Increasing the interactions between atoms (i.e. increasing the scattering length) tends to reduce 
the sound velocity, as can be seen in Fig. HO This behaviour is expected on the grounds that 
interaction strength increases the density of states in general, which increases the compressibility 
of the system, which in turn reduces the speed of sound. Note that in contrast, for a Bose-Einstein 
condensate with a positive scattering length the sound velocity actually increases as interactions 
become stronger. 

V. NON-CUBIC LATTICES AND DIMENSIONAL CROSSOVER 

We now study density response and the speed of sound in non-cubic lattices. By varying the 
lattice heights s x , s y and s z , i.e. by changing the different laser intensities, the cubic symmetry of 
the lattice can be broken. For concreteness, we choose s x = s y . With this simplification, there are 
basically two distinct cases, namely s x > s z and s x < s z . In the former case the atoms are more 
free to move in the z-direction than in the xy-plane. This means that the lattice resembles a square 
lattice arrangement of coupled one-dimensional tubes. The latter corresponds to a situation where 
movement is more restricted in the z-direction. In this case the lattice can thought to be a pile of 
connected two-dimensional planes ("pancakes"). 

In either case, the energy gap A, does not show any qualitative change when the symmetry 
decreases: A merely grows as either s x (and s y ) or s z is increased, since tighter confinement 
effectively increases interaction, see Fig. However, there is a clear signal of a dimensional 
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FIG. 7: The energy gap as a function of the filling fraction with s x = s y < s z , corresponding to quasi-two- 
dimensional pancakes along the xy-plane. The different curves are for different ratios of s x /s z : (from top to 
bottom) 0.5, 0.62, 0.71, 0.83, and 1.0 (cubic case). There is no qualitative difference in the different curves. 
Here the scattering length a is all the time — lOOOao 
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FIG. 8: The speed of sound as a function of the lattice filling fraction in non-cubic lattices with s x = s y < s z , 
corresponding to quasi-two-dimensional pancakes along the iry-plane. The left figure shows the speed in the 
plane, i.e. parallel to the pancakes, and the one on the right the speed orthogonal to the them, i.e. along 
the z-axis. The different curves are for different ratios of s x /s z : (from top to bottom) 1.0 (cubic case), 0.83, 
0.71, 0.62, and 0.5. Here the scattering length a is all the time — 1000ao. 

crossover in the speed of sound. In the two-dimensional limit, i.e. when the tunneling in the 
z-direction is suppressed (pancakes), the local minimum of the speed of sound vanishes, as can be 
seen from Fig. |H1 Also it can be seen that the speed of sound is significantly larger in the direction 
parallel to the planes than orthogonal to them. 

In the opposite case of stronger tunneling in the z-direction than in the xy-plane (ID tubes), 
the local minimum in the speed of sound also disappears. However, as is apparent from Fig. El 
in this one-dimensional limit the speed of sound is again qualitatively different from both the one 
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FIG. 9: The speed of sound as a function of the lattice filling fraction in non-cubic lattices with s x = s y > s z , 
corresponding to quasi-one-dimcnsional tubes along the z-axis. The figure on the left shows the speed in 
the z-direction, i.e. parallel to the tubes, and the right one the speed orthogonal to the them, i.e. in the 
xy-plane. The different curves are for different ratios of s x /s z : (from top to bottom) 1.0 (cubic case), 1.2, 
1.4, 1.6, 1.8, and 2.0. Here the background scattering length a is all the time — 1000a . 

in the three-dimensional cubic system and the one in the non-cubic system in the two-dimensional 
limit. 



In this paper, we have studied the velocity of sound in cubic and non-cubic three-dimensional op- 
tical lattices. We have investigated how the van Hove singularity of the free Fermi gas is smoothened 
by interactions and eventually vanishes when interactions are strong enough. In non-cubic lattices 
we have shown that although the energy gap has a simple behaviour as a function of symmetry, 
the speed of sound shows qualitatively different behaviour over the crossover, and provides a clear 
experimentally observable signature of a dimensional crossover both in the ID and 2D limits. 
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